Measuring the Confidence of the k-Nearest Neighbors Algorithm developed in Python for Visualizing the 3D Stratigraphic Architecture of the Llobregat River Delta¶

In this notebook, we draw maps that measure, in some way, how solid are the predictions used for the horizontal sections in https://www.mdpi.com/2077-1312/10/7/986. For this purpose we take a grid of points in the Llobregat River Delta and we assign to each point a confidence metric satisfying some desirable properties.

First, we configure the directory where the data can be found and also the directory where the figures will be saved.

In [1]:
DATADIR='data/' # Directory with the data
FIGURESDIR='figures/' # Figures produced

Now, we import the packages and the functions that we are going to use. The customer functions are dfined in the file functions.py.

In [2]:
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
import os
In [3]:
from sklearn.neighbors import NearestNeighbors

from sklearn.neighbors._base import _get_weights
In [4]:
import plotly.offline as go_offline
import plotly.graph_objects as go
In [5]:
import functions
from functions import *

The metric we propose uses the weight of the nearest neighbor and divies it by the sum of the weights of the four nearest neighbors such that each one belongs to a different class, this includes the nearest neighbor used to infer the granulometry, so the confidence degree lies always in the interval $[0.25,1]$.

A bonus of confidence is applied when there are several nearest neighbors of the class basement closer than any other non-basement neighbor.

In [6]:
# Global variables needed in the following function
# They are NearestNeighbors objects from sklearn that will be fit before calling the function
knn_sand = NearestNeighbors(n_neighbors=1)
knn_grav = NearestNeighbors(n_neighbors=1) 
knn_clay = NearestNeighbors(n_neighbors=1)
knn_base = NearestNeighbors(n_neighbors=4)

The function confidence(X) will do the job.

In [7]:
def confidence(X):       # The variable `X` is a row of points in the grid.
    
    # For each class, we take the nearest neighbor of each point in X.
    neight_grav=knn_grav.kneighbors(X)[0].reshape((len(X))) 
    neight_sand=knn_sand.kneighbors(X)[0].reshape((len(X)))
    neight_clay=knn_clay.kneighbors(X)[0].reshape((len(X)))
    
    # For the basement, we take up to 4 nearest neighbors to apply the confidence bonus
    m = min(knn_base.n_samples_fit_, 4)
    neight_base=knn_base.kneighbors(X, m)[0].reshape((len(X),m))
        
    conf=np.zeros(len(X))    # Array that will store the confidence degree of each point
        
    for i in range(len(X)): # Here, we conpute the confidence for each point
        bas=sorted(list(neight_base[i]))
        no_bas=[neight_grav[i],neight_sand[i],neight_clay[i]]
        m_no_bas=min(no_bas)
        w=_get_weights(np.array([no_bas+[bas[0]]]), 'distance')[0]
        conf[i]=np.max(w)/np.sum(w)
        
        # Bonus of confidence when several basement points are closer than the nearest non-basement neighbor.
        s=1
        while s<m and bas[s]<m_no_bas:
            w=_get_weights(np.array([no_bas+[bas[s]]]), 'distance')[0]
            conf[i]=conf[i]+(1-conf[i])*(w[3]-max(w[:3]))/np.sum(w)
            s+=1
        
    return 100*conf # We take the confidence to the range 25-100%

The function conf(data,poly,height,axi,grid), were the variables meaning are:

  • data: data points.
  • poly: contour polygon.
  • height: height at which the confidence map will be obtained.
  • axi: axis limits.
  • grid: grid size (number of points in a row/column).

computes the confidence map at a given height. The output of this function is a list [xx,yy,C] where:

  • xx is the X UTM coordinate of a point in the grid.
  • yyis the Y UTM coordinate of a point in the grid.
  • C is a matrix where C[i,j] is the confidence estimated for the point (i,j).
In [8]:
def conf(data,poly,height,axi,grid):
    xx, yy = np.meshgrid(np.linspace(axi[0],axi[1],grid),np.linspace(axi[2],axi[3],grid)) # We create the grid
    C=np.zeros(xx.shape,dtype = float)  # We initalize the matrix
    outside=[]                      # Points outside the contour
    for i in range(xx.shape[0]):
        for j in range(xx.shape[1]):
            if poly.distance(geometry.Point(xx[i,j],yy[i,j]))>1:
                outside.append((i,j))
    
    df=data[data.Cota==height] # We take the data points at the given height
        
    # We split the data into the four categories
    gravels=df[df.Clase=='gravillas y gravas']
    sands=df[df.Clase=='arenas']
    clays=df[df.Clase=='arcillas y limos']
    basements=df[df.Valor=='S']
    
    # We fit the global NearestNeighbors objects
    knn_grav.fit(list(zip(gravels['UTM_X'],gravels['UTM_Y'])))
    knn_sand.fit(list(zip(sands['UTM_X'],sands['UTM_Y'])))
    knn_clay.fit(list(zip(clays['UTM_X'],clays['UTM_Y'])))
    knn_base.fit(list(zip(basements['UTM_X'],basements['UTM_Y'])))
    
    # We call the confidence function
    for i in range(xx.shape[0]):
        d=list(zip(xx[i],yy[i]))
        C[i]=confidence(d)

    for (i,j) in outside: # We set the confidence to nan outside the contour
        C[i,j]=np.nan
        
    return [xx,yy,C]

We decide to create a confident map every 5 meters depth.

In [9]:
Heights=list(range(0,-101,-5)) # Heights where the confidence maps will be obtained

We also set the plot frame.

In [10]:
FRAME=500 # Plot frame

We read the data.

In [11]:
boreholes_data=pd.read_excel(DATADIR+'hsd new basements.xls') # We read the data from the xls sheet
boreholes_data=boreholes_data.drop(columns=['Codi'])  # We do not need this column

ax=[boreholes_data.UTM_X.min()-FRAME,
    boreholes_data.UTM_X.max()+FRAME,
    boreholes_data.UTM_Y.min()-FRAME,
    boreholes_data.UTM_Y.max()+FRAME] # Axis limits (with frame)
minx_rounded = 1000 * round(ax[0]/1000)  # We round the X coordinate limits
maxx_rounded = 1000 * round(ax[1]/1000) 
In [12]:
contour=pd.read_csv(DATADIR+'deltacontour.csv')
contour=contour.drop(columns=['Cota'])          
contour_poly=geometry.Polygon(zip(contour['UTM_X'],contour['UTM_Y']))  # Delta polygon
xpol,ypol = contour_poly.exterior.xy # Delta contour coordinates
In [13]:
xx,yy,C=conf(boreholes_data,contour_poly,-100,ax,200) # Confidence map at 100m b.s.l. (below sea level)

The function coordinates creates a list with three lists (UTM_X, UTM_Y and Height coordinates).

In [14]:
xyzcontour=coordinates(DATADIR+'deltacontour.csv',[0,1,2]) 
In [15]:
bound=bounds(xyzcontour) # bounds for the figures

The funtion conf_fig(height) creates a figure with the confidence index at a given heigh.

In [16]:
def conf_fig(height):
    # Confidence ranges
    cmap = plt.cm.rainbow
    cmaplist = [cmap(i) for i in range(cmap.N)]
    cmap = mpl.colors.LinearSegmentedColormap.from_list('Custom cmap', cmaplist, cmap.N)
    bounds=np.array([25, 30, 35, 40, 50, 60, 70, 80, 90, 95, 100])
    norm1=mpl.colors.BoundaryNorm(bounds, cmap.N)
    
    # figure data
    xx,yy,C=conf(boreholes_data,contour_poly,height,ax,200)
    
    # create de figure
    fig, ax1 = plt.subplots()
    
    ax1.imshow(C, cmap=cmap, origin='lower', norm=norm1, alpha=0.6,
           extent=[xx.min(), xx.max(), yy.min(), yy.max()])
    
    sm = plt.cm.ScalarMappable(norm=norm1, cmap=cmap)
    sm.set_array([])
    
    fig.colorbar(sm,spacing='proportional').set_label(' %', rotation=0)
    
    ax1.plot(xpol, ypol, alpha=0.6, color='black', linewidth=1.5)    

    ax1.set_xlabel('UTM_X')
    ax1.set_ylabel('UTM_Y')

    ax1.set_xticks(np.arange(minx_rounded, maxx_rounded, step=3000))
    
    ax1.set_title("Height "+ str(height) +' m')
    ax1.axis(ax)

    return ax1
    

For example, the confidence map at height -10 is

In [17]:
conf_fig(-10)
Out[17]:
<AxesSubplot:title={'center':'Height -10 m'}, xlabel='UTM_X', ylabel='UTM_Y'>

The function knn_fig(h) uses the funtion knn we defined in the paper https://www.mdpi.com/2077-1312/10/7/986, which can be found explicitly in function.py, to draw a figure with the knn precictions at height h.

In [18]:
def knn_fig(h):
    # figure settings and bounds: we set a color for each category
    SCATTER=True
    OPACITY=0.6
    cm = mpl.colors.ListedColormap(['cyan','black','yellow','black','gray','black','saddlebrown'])
    cm = cm(np.arange(cm.N))
    cm[:,-1]=np.array([OPACITY,0,OPACITY,0,OPACITY,0,OPACITY])
    cm = mpl.colors.ListedColormap(cm)
    epsilon=10**(-15)
    levels = [-0.5,0.5,1-epsilon,1+epsilon,1.5,2+epsilon,2.5,3.5]
    norm=mpl.colors.BoundaryNorm(levels, 7)
    
    # figure data
    layer=layer_function(boreholes_data,h)
    xyC=knn(boreholes_data,contour_poly,h,ax,200)
    
    # create de figure
    fig2, ax2 = plt.subplots()
    
    ax2.contourf(xyC[0],xyC[1],xyC[2], levels, cmap=cm, norm=norm, antialiased = True)
    
    if SCATTER:
        ax2.scatter(layer[1]['UTM_X'], layer[1]['UTM_Y'], c='blue', s=5, alpha=0.8, label='gravels')
        ax2.scatter(layer[2]['UTM_X'], layer[2]['UTM_Y'], c='orange', s=5, alpha=0.7, label='sands')
        ax2.scatter(layer[3]['UTM_X'], layer[3]['UTM_Y'], c='black', s=5, alpha=0.6, label='clays')
        ax2.scatter(layer[4]['UTM_X'], layer[4]['UTM_Y'], c='brown', s=5, alpha=0.6, label='basement')

    ax2.plot(xpol, ypol, alpha=0.6, color='black', linewidth=1.5) 
    ax2.set_xlabel('UTM_X')
    ax2.set_ylabel('UTM_Y')

    ax2.set_xticks(np.arange(minx_rounded, maxx_rounded, step=3000))

    ax2.set_title("Height "+str(h)+' m')
    ax2.axis(ax)
    return ax2
    

For example the knn prediction figure at height -10 is

In [19]:
knn_fig(-10)
Out[19]:
<AxesSubplot:title={'center':'Height -10 m'}, xlabel='UTM_X', ylabel='UTM_Y'>

We define the funtion cyk_fig(h) to create a figure with the knn predictions and the confidence maps together at a given height h.

In [20]:
def cyk_fig(h):
    
    cmap = plt.cm.rainbow
    cmaplist = [cmap(i) for i in range(cmap.N)]
    cmap = mpl.colors.LinearSegmentedColormap.from_list('Custom cmap', cmaplist, cmap.N)
    bounds=np.array([25, 30, 35, 40, 50, 60, 70, 80, 90, 95, 100])
    norm1=mpl.colors.BoundaryNorm(bounds, cmap.N)
    
    SCATTER=True
    OPACITY=0.6
    cm = mpl.colors.ListedColormap(['cyan','black','yellow','black','gray','black','saddlebrown'])
    cm = cm(np.arange(cm.N))
    cm[:,-1]=np.array([OPACITY,0,OPACITY,0,OPACITY,0,OPACITY])
    cm = mpl.colors.ListedColormap(cm)
    epsilon=10**(-15)
    levels = [-0.5,0.5,1-epsilon,1+epsilon,1.5,2+epsilon,2.5,3.5]
    norm2=mpl.colors.BoundaryNorm(levels, 7)
    
    # figure data
    xx,yy,C=conf(boreholes_data,contour_poly,h,ax,300)
    
    layer=layer_function(boreholes_data,h)
    xyC=knn(boreholes_data,contour_poly,h,ax,300)
    
    # figure creation 
    fig = plt.figure(figsize=(14, 4))
    
    # Adds subplot on position 2
    ax1 = fig.add_subplot(122)
    # Adds subplot on position 1
    ax2 = fig.add_subplot(121)
    
    ax1.imshow(C, cmap=cmap, origin='lower', norm=norm1, alpha=0.6, aspect="auto",
           extent=[xx.min(), xx.max(), yy.min(), yy.max()])
    
    
    ax1.plot(xpol, ypol, alpha=0.6, color='black', linewidth=1.5)    

    ax1.set_xlabel('UTM_X')
    #ax1.set_ylabel('UTM_Y')

    ax1.set_xticks(np.arange(minx_rounded, maxx_rounded, step=3000))
    
    ax1.set_title("Height "+ str(h) +' m')
    ax1.axis(ax)
    ax1.set_yticklabels([])    # Removes Y axis labels
    
    ax2.contourf(xyC[0],xyC[1],xyC[2], levels, cmap=cm, norm=norm2, antialiased = True)
    
    if SCATTER:
        ax2.scatter(layer[1]['UTM_X'], layer[1]['UTM_Y'], c='blue', s=5, alpha=0.8, label='gravels')
        ax2.scatter(layer[2]['UTM_X'], layer[2]['UTM_Y'], c='orange', s=5, alpha=0.7, label='sands')
        ax2.scatter(layer[3]['UTM_X'], layer[3]['UTM_Y'], c='black', s=5, alpha=0.6, label='clays')
        ax2.scatter(layer[4]['UTM_X'], layer[4]['UTM_Y'], c='brown', s=5, alpha=0.6, label='basement')

    ax2.plot(xpol, ypol, alpha=0.6, color='black', linewidth=1.5) 
    ax2.set_xlabel('UTM_X')
    ax2.set_ylabel('UTM_Y')

    ax2.set_xticks(np.arange(minx_rounded, maxx_rounded, step=3000))

    ax2.set_title("Height "+str(h)+' m')
    ax2.axis(ax)
    
    sm = plt.cm.ScalarMappable(norm=norm1, cmap=cmap)
    sm.set_array([])
    
    fig.colorbar(sm,spacing='proportional', ax=[ax2,ax1]).set_label(' %', rotation=0)
    
    return (fig,ax1,ax2)
    

For example at height -10 we get

In [21]:
cyk_fig(-10)
Out[21]:
(<Figure size 1008x288 with 3 Axes>,
 <AxesSubplot:title={'center':'Height -10 m'}, xlabel='UTM_X'>,
 <AxesSubplot:title={'center':'Height -10 m'}, xlabel='UTM_X', ylabel='UTM_Y'>)

We create a list with all figures every 5m.

In [22]:
figures=[cyk_fig(h) for h in Heights]
/var/folders/fd/bd2ntp3n1lq8rv9w3gj8d8gw0000gn/T/ipykernel_68207/739094612.py:26: RuntimeWarning:

More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).

And we save the figures in the figures directory.

In [23]:
for i in range(21):
    
    filename='height'+str(-5*i)+'.png'
    
    filename=FIGURESDIR+filename
    os.makedirs(os.path.dirname(filename), exist_ok=True)
    figures[i][0].savefig(filename, transparent=False, dpi=164, bbox_inches='tight')

3D figures¶

Basement¶

The confidence for basement. We use a $200\times200$ grid.

In [24]:
# knn algorith predictions for a square grid of width 200

data_split=[[],[],[],[]]
for h in Heights:
    cf=conf(boreholes_data,contour_poly,h,ax,200)
    cc=knn(boreholes_data,contour_poly,h,ax,200)
    for c in range(4):
        x=[cc[0][i,j] for i in range(cc[0].shape[0]) for j in range(cc[0].shape[1]) if cc[2][i,j]==c]
        y=[cc[1][i,j] for i in range(cc[1].shape[0]) for j in range(cc[1].shape[1]) if cc[2][i,j]==c]
        z=[h for i in range(cc[0].shape[0]) for j in range(cc[0].shape[1]) if cc[2][i,j]==c]
        ccf=[cf[2][i,j] for i in range(cc[1].shape[0]) for j in range(cc[1].shape[1]) if cc[2][i,j]==c]
        data_split[c].append([x,y,z,ccf])

gravels_knn200=data_split[0] # knn predictions for gravels
sands_knn200=data_split[1]  # knn predictions for sands
clays_knn200=data_split[2]  # knn predictions for clays
basement_knn200=data_split[3] # knn predictions for basement

We group the confidence and the coordinates, getting four different lists.

In [25]:
x_bas200=np.concatenate([basement_knn200[i][0] for i in range(21)]) # X UTM coordinates
y_bas200=np.concatenate([basement_knn200[i][1] for i in range(21)]) # Y UTM coordinates
z_bas200=np.concatenate([basement_knn200[i][2] for i in range(21)]) # Z coordinates
conf_bas200=np.concatenate([basement_knn200[i][3] for i in range(21)]) # Confidence
In [26]:
xyzc_bas200=list(zip(x_bas200,y_bas200,z_bas200,conf_bas200)) # We zip the four lists
In [27]:
xyz_bas200=[x[:3] for x in xyzc_bas200] # Here, we take only the XYZ coordinates

We can work out the average confidence of the basement point for a map. Also the maximum and minimum heights where basement can be found and the medium height of the basement points.

In [28]:
bas_conf_mean=np.mean(conf_bas200)
bas_height_max=np.max(z_bas200)
bas_height_min=np.min(z_bas200)
bas_height_mean=np.mean(z_bas200)
In [29]:
(bas_conf_mean,bas_height_max,bas_height_min,bas_height_mean)
Out[29]:
(77.98397866313435, 0, -100, -70.41003196387203)

The basement surface¶

Now, we set the grid size to 50, it is big enough and allows us to perform the computation within a reasonable time.

In [30]:
basementknn50max=pd.read_csv(DATADIR+'basementknn50max.csv')
In [31]:
basementknn50max
Out[31]:
UTM_X UTM_Y Cota
0 426262 4580059 0
1 426640 4580059 0
2 426640 4580345 0
3 427019 4580345 0
4 426640 4580632 0
... ... ... ...
11731 419445 4582063 -35
11732 419824 4582063 -35
11733 418688 4582350 -35
11734 419067 4582350 -35
11735 419445 4582350 -35

11736 rows × 3 columns

In [32]:
basement50_coordinates=coordinates(DATADIR+'basementknn50max.csv',[0,1,2])
In [33]:
# We change the limits taking into account the data in basement_coordinates
basement50_bounds=[min(basement50_coordinates[0]),max(basement50_coordinates[0]),
                 min(basement50_coordinates[1]),max(basement50_coordinates[1]),
                 max(basement50_coordinates[0])-min(basement50_coordinates[0]),
                 max(basement50_coordinates[1])-min(basement50_coordinates[1])]
new_bounds50=bounds_join(bound,basement50_bounds)

We now interpolate

In [34]:
interpolation_knnbasement50=interpolation(basement50_coordinates,50,new_bounds50)

and reduce the data to the points inside the contour:

In [35]:
interpolation_knnbasement50cutted=cutting(interpolation_knnbasement50,contour_poly,500)

We shape the data as we need

In [36]:
xyz_interpolation_knnbasementcutted50=[[interpolation_knnbasement50cutted[1][i],
                                      interpolation_knnbasement50cutted[2][i],
                                      interpolation_knnbasement50cutted[0][i]] for i in range(len(interpolation_knnbasement50cutted[0]))]

Data for first reached basement¶

We set the grid side to 200.

The results from the following four cells are stored in a csv file to save computation time.

Running them is not mandatory if the file basementknn200confidence.csv is present in the data folder.

The reader may just load the data as shown further.

In [ ]:
# Data for the points where the basements is first reached

basementknn200max=pd.read_csv(DATADIR+'basementknn200max.csv')

#  x,y,z list for the data in basementknn200max

basement200max_coordinates=coordinates(DATADIR+'basementknn200max.csv',[0,1,2])

# We change the limits taking into account the data in basement_coordinates
basement200_bounds=[min(basement200max_coordinates[0]),max(basement200max_coordinates[0]),
                 min(basement200max_coordinates[1]),max(basement200max_coordinates[1]),
                 max(basement200max_coordinates[0])-min(basement200max_coordinates[0]),
                 max(basement200max_coordinates[1])-min(basement200max_coordinates[1])]
new_bounds200=bounds_join(bound,basement200_bounds)
In [ ]:
xyz_basement200max=list(zip(basement200max_coordinates[0],
                            basement200max_coordinates[1],
                            basement200max_coordinates[2]))
In [ ]:
xyzc_basement200max=[x for x in xyzc_bas200 if ((int(x[:3][0]),int(x[:3][1]),int(x[:3][2])) in  xyz_basement200max) ]
In [ ]:
# We save the reduced data to a csv file to save computation time
df=pd.DataFrame(xyzc_basement200max,columns=["UTM_X","UTM_Y","Cota","Confidence"])
df.to_csv(DATADIR+'basementknn200confidence.csv',index=False)

Here, the csv file is loaded.

In [37]:
bknn200cmax=pd.read_csv(DATADIR+'basementknn200confidence.csv')
In [38]:
bknn200cmax
Out[38]:
UTM_X UTM_Y Cota Confidence
0 426062.472362 4.579741e+06 0 28.593293
1 426155.718593 4.579741e+06 0 30.005940
2 425969.226131 4.579812e+06 0 28.314688
3 426062.472362 4.579812e+06 0 31.066904
4 426155.718593 4.579812e+06 0 33.898017
... ... ... ... ...
12587 417017.587940 4.569730e+06 -100 35.940012
12588 417110.834171 4.569730e+06 -100 35.901846
12589 417204.080402 4.569730e+06 -100 35.360105
12590 416924.341709 4.569801e+06 -100 33.214541
12591 417017.587940 4.569801e+06 -100 33.599572

12592 rows × 4 columns

In [39]:
xyzc_b200max=coordinatesc(DATADIR+'basementknn200confidence.csv',[0,1,2,3])

Now we are ready to draw a figure with the basement surface and the confidence of the prediction of each point where basement is first predicted.

In [40]:
fig=go.Figure()

fig.add_trace(go.Scatter3d(x=xyzc_b200max[0], y=xyzc_b200max[1], z=xyzc_b200max[2],
            mode ='markers',
            name='knn_basament',
            marker = dict(size = 2,
                          color = xyzc_b200max[3],
                          colorscale='rainbow',
                          opacity = 1,
                          symbol='circle',
                         showscale=True)
                        ))

fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], 
                           mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )))

fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
                         x=interpolation_knnbasement50cutted[1],
                         y=interpolation_knnbasement50cutted[2], 
                        opacity = 0.5,
                        #inensity=interpolation_knnbasementcutted[0],
                        colorscale='oryel', 
                        name='basement surface',
                        showscale=False
                        ))

fig.update_layout(legend=dict(
    yanchor="top",
    y=0.99,
    xanchor="left",
    x=0.01))

fig.update_layout( title="3D Basement surface Llobregat Delta and confidence index, Z scale is x 50.",
    scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
                             xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
                             yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
                            ))


fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_Basement_Confidence.html',validate=True, auto_open=False)
Out[40]:
'figures/3D_Basement_Confidence.html'

3D figure with the confidence of the basement at each heigh (without the basement surface).

In [41]:
fig=go.Figure()

fig.add_trace(go.Scatter3d(x=x_bas200, y=y_bas200, z=z_bas200,
            mode ='markers',
            name='knn_basament',
            marker = dict(size = 2,
                          color = conf_bas200,
                          colorscale='rainbow',
                          opacity = 1,
                          symbol='circle',
                         showscale=True)
                        ))

fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], 
                           mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )))

#fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
 #                        x=interpolation_knnbasement50cutted[1],
  #                       y=interpolation_knnbasement50cutted[2], 
   #                     opacity = 0.5,
    #                    #inensity=interpolation_knnbasementcutted[0],
     #                   colorscale='oryel', 
      #                  name='basement surface',
       #                 showscale=False
        #                ))

fig.update_layout(legend=dict(
    yanchor="top",
    y=0.99,
    xanchor="left",
    x=0.01))

fig.update_layout( title="3D Basement surface Llobregat Delta and confidence index all height, Z scale is x 50.",
    scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
                             xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
                             yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
                            ))


fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'Basement_Confidence_Layers.html',validate=True, auto_open=False)
Out[41]:
'figures/Basement_Confidence_Layers.html'

For gravels and sands¶

Next, we will obtain two figures. The first one shows the gravel predictions colored by confidence. The second one shows the same for the sand predictions.

Confidence list.

In [42]:
conf_gr=np.concatenate([gravels_knn200[i][3] for i in range(21)])
conf_sn=np.concatenate([sands_knn200[i][3] for i in range(21)])

Averages.

In [43]:
gr200_mean=np.mean(conf_gr)
sn200_mean=np.mean(conf_sn)
In [44]:
gr200_mean # Average confidence for gravel predictions
Out[44]:
49.33703116232511
In [45]:
sn200_mean # Average confidence for sand predictions
Out[45]:
50.18378710802425
In [46]:
x_gr=np.concatenate([gravels_knn200[i][0] for i in range(21)]) # X coordiates of gravel predictions
y_gr=np.concatenate([gravels_knn200[i][1] for i in range(21)]) # Y coordiates of gravel predictions
z_gr=np.concatenate([gravels_knn200[i][2] for i in range(21)]) # Z coordiates of gravel predictions

The following code cell generate the figure with the confidence of the gravel predictions.

In [47]:
fig=go.Figure()
x=x_gr
y=y_gr
z=z_gr
c=conf_gr

fig.add_trace(go.Scatter3d(x=x, y=y, z=z,
            mode ='markers',
            name='knn_gravels',
            marker = dict(size = 2,
                          color = c,
                          colorscale='rainbow',
                          opacity = 1,
                          symbol='circle',
                         showscale=True)
                        ))

fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], 
                           mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )))

fig.update_layout(legend=dict(
    yanchor="top",
    y=0.99,
    xanchor="left",
    x=0.01))

fig.update_layout( title="3D confidence map for gravels index, Z scale is x 50.",
    scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
                             xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
                             yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
                            ))


fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'Gravel_Lithosomes_Confidence_Layers.html',validate=True, auto_open=False)
Out[47]:
'figures/Gravel_Lithosomes_Confidence_Layers.html'

Now, we do the same for sand predictions.

In [48]:
x_sn=np.concatenate([sands_knn200[i][0] for i in range(21)])
y_sn=np.concatenate([sands_knn200[i][1] for i in range(21)])
z_sn=np.concatenate([sands_knn200[i][2] for i in range(21)])

fig=go.Figure()
x=x_sn
y=y_sn
z=z_sn
c=conf_sn

fig.add_trace(go.Scatter3d(x=x, y=y, z=z,
            mode ='markers',
            name='knn_sands',
            marker = dict(size = 2,
                          color = c,
                          colorscale='rainbow',
                          opacity = 1,
                          symbol='circle',
                         showscale=True)
                        ))

fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], 
                           mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )))

fig.update_layout(legend=dict(
    yanchor="top",
    y=0.99,
    xanchor="left",
    x=0.01))

fig.update_layout( title="3D confidence map for sands index, Z scale is x 50.",
    scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
                             xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
                             yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
                            ))


fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'Sand_Lithosomes_Confidence_Layers.html',validate=True, auto_open=False)
Out[48]:
'figures/Sand_Lithosomes_Confidence_Layers.html'

Lithosomes¶

We shall plot gravel and sand lithosomes like in the previous work https://www.mdpi.com/2077-1312/10/7/986. Then, we will study the confience of the predictions in each lithosome. We shall display tables with the average confidence and draw the lithosomes colored by average confidence.

In [49]:
# knn algorith predictions for a grid of size 50

data_split=[[],[],[],[]]
for h in Heights:
    cf=conf(boreholes_data,contour_poly,h,ax,50)
    cc=knn(boreholes_data,contour_poly,h,ax,50)
    for c in range(4):
        x=[cc[0][i,j] for i in range(cc[0].shape[0]) for j in range(cc[0].shape[1]) if cc[2][i,j]==c]
        y=[cc[1][i,j] for i in range(cc[1].shape[0]) for j in range(cc[1].shape[1]) if cc[2][i,j]==c]
        z=[h for i in range(cc[0].shape[0]) for j in range(cc[0].shape[1]) if cc[2][i,j]==c]
        ccf=[cf[2][i,j] for i in range(cc[1].shape[0]) for j in range(cc[1].shape[1]) if cc[2][i,j]==c]
        data_split[c].append([x,y,z,ccf])

gr_knn50=data_split[0] # knn predictions for gravels
sn_knn50=data_split[1]  # knn predictions for sands
In [50]:
# We split X, Y, Z and confidence. For gravels
x_gr50=np.concatenate([gr_knn50[i][0] for i in range(21)]) 
y_gr50=np.concatenate([gr_knn50[i][1] for i in range(21)]) 
z_gr50=np.concatenate([gr_knn50[i][2] for i in range(21)])
c_gr50=np.concatenate([gr_knn50[i][3] for i in range(21)])
In [51]:
xyzc_gr50=list(zip(x_gr50,y_gr50,z_gr50,c_gr50))
In [52]:
# Now for sands
x_sn50=np.concatenate([sn_knn50[i][0] for i in range(21)])
y_sn50=np.concatenate([sn_knn50[i][1] for i in range(21)])
z_sn50=np.concatenate([sn_knn50[i][2] for i in range(21)])
c_sn50=np.concatenate([sn_knn50[i][3] for i in range(21)])
In [53]:
xyzc_sn50=list(zip(x_sn50,y_sn50,z_sn50,c_sn50))
In [54]:
# We compute the average confidence
gr50_mean=np.mean(c_gr50) # for gravel predictions
sn50_mean=np.mean(c_sn50) # for sand predictions
In [55]:
gr50_mean,sn50_mean
Out[55]:
(49.57830307839629, 50.20226163638911)

We compare with those obtained with grid size=200. There is not much difference.

In [56]:
gr200_mean,sn200_mean
Out[56]:
(49.33703116232511, 50.18378710802425)

We compute the lithosomes as in https://www.mdpi.com/2077-1312/10/7/986.

In [57]:
# we take the x ,y and z coordinates for the knn  predictions

gravels_knn50=[xyc(boreholes_data,contour_poly,ax,0,h,50) for h in Heights]
sands_knn50=[xyc(boreholes_data,contour_poly,ax,1,h,50) for h in Heights]
In [58]:
# Points for gravels.
# p1=(414522,4568893,-115) we are taking max depth = -100, so this lithosome is removed 
p2=(418355, 4570004,-80)
p3=(423682,4573012,-60)
p4=(424747,4574333,-65)
p5=(420960, 4571183,-65)
p6=(427019,4578055,-50)
p7=(422096,4577769,-30)
p8=(414901,4570897,-30)
p9=(417173,4571183,-30)
p10=(427777,4579773,-30)
p11=(430427,4579773,-25)
p12=(418309,4573474,-20)
p13=(417931,4582064,-10)
p14=(425126,4577769,0)

gravels_selection1=[#p1,
                    p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14]
gravels_selection2=[[x[0] for x in gravels_selection1],
                    [x[1] for x in gravels_selection1],
                    [x[2] for x in gravels_selection1]]

# Points for sands.
q1=(423611,4570611,-100)
q2=(415658,4568606,-90)
q3=(422096,4570038,-85)
q4=(425883,4572042,-80)
q5=(427398,4576623,-70)
q6=(423990,4575192,-55)
q7=(414144,4568893,-15)
q8=(418688,4570611,-10)
q9=(419445,4574903,-5)
q10=(428155,4576910,-5)
q11=(420203,4572901,0)
q12=(423232,4572329,0)
q13=(422475,4576051,0)
q14=(428155,4574333,0)
q15=(425883,4578628,0)
q16=(427398,4578055,0)
q17=(430806,4580059,0)

sands_selection1=[q1,q2,q3,q4,q5,q6,q7,q8,q9,q10,q11,q12,q13,q14,q15,q16,q17]
sands_selection2=[[x[0] for x in sands_selection1],
                    [x[1] for x in sands_selection1],
                    [x[2] for x in sands_selection1]]
In [59]:
# lithosomes for gravels

xygr_knn50=[list(zip(x[0],x[1])) for x in gravels_knn50]

#gravels_lit1=lithosome(p1,xygr_knn50,400,Heights,'grlit1')
#gravels_lit1_b=[x for x in gravels_lit1[0] if x[2]< -35]
#h1_b = ConvexHull(gravels_lit1_b)
#gravels_lit1b=[h1_b.points,h1_b.vertices,h1_b.simplices,'grlit1']

gravels_lit2=lithosome(p2,xygr_knn50,400,Heights,'grlit1')
gravels_lit3=lithosome(p3,xygr_knn50,400,Heights,'grlit2')
gravels_lit4=lithosome(p4,xygr_knn50,400,Heights,'grlit3')
gravels_lit5=lithosome(p5,xygr_knn50,400,Heights,'grlit4')
gravels_lit6=lithosome(p6,xygr_knn50,400,Heights,'grlit5')
gravels_lit7=lithosome(p7,xygr_knn50,400,Heights,'grlit6')
gravels_lit8=lithosome(p8,xygr_knn50,400,Heights,'grlit7')
gravels_lit9=lithosome(p9,xygr_knn50,400,Heights,'grlit8')
gravels_lit10=lithosome(p10,xygr_knn50,400,Heights,'grlit9')
gravels_lit11=lithosome(p11,xygr_knn50,500,Heights,'grlit10')
gravels_lit12=lithosome(p12,xygr_knn50,500,Heights,'grlit11')
gravels_lit13=lithosome(p13,xygr_knn50,400,Heights,'grlit12')
gravels_lit14=lithosome(p14,xygr_knn50,400,Heights,'grlit13')


# List of gravels lithosomes 

list_grlithosomes=[#gravels_lit1b,
                   gravels_lit2,gravels_lit3,gravels_lit4,
                  gravels_lit5,gravels_lit6,gravels_lit7,gravels_lit8,
                  gravels_lit9,gravels_lit10,gravels_lit11,gravels_lit12,
                  gravels_lit13,gravels_lit14]
In [60]:
# lithosome for sands

xysnd_knn=[list(zip(x[0],x[1])) for x in sands_knn50]

snd_lit1=lithosome(q1,xysnd_knn,400,Heights,'sndlit1')
snd_lit1_b=[x for x in snd_lit1[0] if x[2]< -89]
ha1_1b = ConvexHull(snd_lit1_b)
snd_lit1b=[ha1_1b.points,ha1_1b.vertices,ha1_1b.simplices,'sndlit1']

snd_lit2=lithosome(q2,xysnd_knn,400,Heights,'sndlit2')
snd_lit2_b=[x for x in snd_lit2[0] if x[2]<-70]
ha2_1b = ConvexHull(snd_lit2_b)
snd_lit2b=[ha2_1b.points,ha2_1b.vertices,ha2_1b.simplices,'sndlit2']

snd_lit3=lithosome(q3,xysnd_knn,400,Heights,'sndlit3')
snd_lit3_b=[x for x in snd_lit3[0] if x[2]< -70 and -90<x[2]]+[ 
 [ 4.20960571e+05,  4.56975163e+06, -9.00000000e+01],
 [ 4.21339265e+05,  4.56946531e+06, -9.00000000e+01],
 [ 4.21339265e+05,  4.56975163e+06, -9.00000000e+01],
 [ 4.21339265e+05,  4.57003796e+06, -9.00000000e+01],
 [ 4.21717959e+05,  4.56975163e+06, -9.00000000e+01],
 [ 4.21717959e+05,  4.57003796e+06, -9.00000000e+01],
 [ 4.21717959e+05,  4.57032429e+06, -9.00000000e+01],
 [ 4.21717959e+05,  4.57061061e+06, -9.00000000e+01],
 [ 4.22096653e+05,  4.56975163e+06, -9.00000000e+01],
 [ 4.22096653e+05,  4.57003796e+06, -9.00000000e+01],
 [ 4.22096653e+05,  4.57032429e+06, -9.00000000e+01],
 [ 4.22096653e+05,  4.57061061e+06, -9.00000000e+01],
 [ 4.22475347e+05,  4.57003796e+06, -9.00000000e+01],
 [ 4.22475347e+05,  4.57032429e+06, -9.00000000e+01],
 [ 4.22475347e+05,  4.57061061e+06, -9.00000000e+01]]
ha3_1b = ConvexHull(snd_lit3_b)
snd_lit3b=[ha3_1b.points,ha3_1b.vertices,ha3_1b.simplices,'sndlit3']

snd_lit4=lithosome(q4,xysnd_knn,400,Heights,'sndlit4')
snd_lit4_b=[x for x in snd_lit4[0] if x[2]< -10 and x[2]>-90]+[
 [ 4.25504898e+05,  4.57118327e+06, -9.00000000e+01],
 [ 4.25504898e+05,  4.57146959e+06, -9.00000000e+01],
 [ 4.25504898e+05,  4.57175592e+06, -9.00000000e+01],
 [ 4.25504898e+05,  4.57204224e+06, -9.00000000e+01],
 [ 4.25504898e+05,  4.57232857e+06, -9.00000000e+01],
 [ 4.25504898e+05,  4.57261490e+06, -9.00000000e+01],
 [ 4.25883592e+05,  4.57232857e+06, -9.00000000e+01],
 [ 4.25883592e+05,  4.57261490e+06, -9.00000000e+01],
 [ 4.25883592e+05,  4.57290122e+06, -9.00000000e+01],
 [ 4.25883592e+05,  4.57318755e+06, -9.00000000e+01],
 [ 4.25883592e+05,  4.57347388e+06, -9.00000000e+01],
 [ 4.26262286e+05,  4.57175592e+06, -9.00000000e+01],
 [ 4.26262286e+05,  4.57204224e+06, -9.00000000e+01],
 [ 4.26262286e+05,  4.57232857e+06, -9.00000000e+01],
 [ 4.26262286e+05,  4.57261490e+06, -9.00000000e+01],
 [ 4.26262286e+05,  4.57290122e+06, -9.00000000e+01],
 [ 4.26262286e+05,  4.57318755e+06, -9.00000000e+01],
 [ 4.2664098e+05,  4.5726149e+06, -9.0000000e+01],
 [ 4.26640980e+05,  4.57290122e+06, -9.00000000e+01]]
ha4_1b = ConvexHull(snd_lit4_b)
snd_lit4b=[ha4_1b.points,ha4_1b.vertices,ha4_1b.simplices,'sndlit4']

snd_lit5=lithosome(q5,xysnd_knn,400,Heights,'sndlit5')
snd_lit5_b=[x for x in snd_lit5[0] if x[2]<-50]
ha5_1b = ConvexHull(snd_lit5_b)
snd_lit5b=[ha5_1b.points,ha5_1b.vertices,ha5_1b.simplices,'sndlit5']

snd_lit6=lithosome(q6,xysnd_knn,400,Heights,'sndlit6')
snd_lit6_b=[x for x in snd_lit6[0] if x[2]<-50]
ha6_1b = ConvexHull(snd_lit6_b)
snd_lit6b=[ha6_1b.points,ha6_1b.vertices,ha6_1b.simplices,'sndlit6']

snd_lit7=lithosome(q7,xysnd_knn,400,Heights,'sndlit7')
snd_lit7_b=[x for x in snd_lit7[0] if x[2]>-30]
ha7_1b = ConvexHull(snd_lit7_b)
snd_lit7b=[ha7_1b.points,ha7_1b.vertices,ha7_1b.simplices,'sndlit8']

snd_lit8=lithosome(q8,xysnd_knn,400,Heights,'sndlit8')
snd_lit8_b=[x for x in snd_lit8[0] if x[2]>-30 and x[2]<-4]
ha8_1b = ConvexHull(snd_lit8_b)
snd_lit8b=[ha8_1b.points,ha8_1b.vertices,ha8_1b.simplices,'sndlit8']

snd_lit9=lithosome(q9,xysnd_knn,400,Heights,'sndlit9')
snd_lit9_b=[x for x in snd_lit9[0] if x[2]>-25]
ha9_1b = ConvexHull(snd_lit9_b)
snd_lit9b=[ha9_1b.points,ha9_1b.vertices,ha9_1b.simplices,'sndlit9']

snd_lit10=lithosome(q10,xysnd_knn,400,Heights,'sndlit10')
snd_lit10_b=[x for x in snd_lit10[0] if x[2]>-35]
ha10_2b = ConvexHull(snd_lit10_b)
snd_lit10b=[ha10_2b.points,ha10_2b.vertices,ha10_2b.simplices,'sndlit10']

snd_lit11=lithosome(q11,xysnd_knn,400,Heights,'sndlit11')
snd_lit11_b=[x for x in snd_lit11[0] if x[2]>-36]
ha11_1b = ConvexHull(snd_lit11_b)
snd_lit11b=[ha11_1b.points,ha11_1b.vertices,ha11_1b.simplices,'sndlit11']

snd_lit12=lithosome(q12,xysnd_knn,400,Heights,'sndlit12')
snd_lit12_b=[x for x in snd_lit12[0] if x[2]>-37 and x[1]<4576053]
ha12_2b = ConvexHull(snd_lit12_b)
snd_lit12b=[ha12_2b.points,ha12_2b.vertices,ha12_2b.simplices,'sndlit12']

snd_lit13=lithosome(q13,xysnd_knn,400,Heights,'sndlit13')
snd_lit13_b=[x for x in snd_lit13[0] if x[2]>-27]
ha13_3b = ConvexHull(snd_lit13_b)
snd_lit13b=[ha13_3b.points,ha13_3b.vertices,ha13_3b.simplices,'sndlit13']

snd_lit14=lithosome(q14,xysnd_knn,400,Heights,'sndlit14')
snd_lit14_b=[x for x in snd_lit14[0] if x[2]>-57]
ha14_4b = ConvexHull(snd_lit14_b)
snd_lit14b=[ha14_4b.points,ha14_4b.vertices,ha14_4b.simplices,'sndlit14']

snd_lit15=lithosome(q15,xysnd_knn,400,Heights,'sndlit15')
snd_lit15_b=[x for x in snd_lit15[0] if  x[1]>4578055 ]
ha15_5b = ConvexHull(snd_lit15_b)
snd_lit15b=[ha15_5b.points,ha15_5b.vertices,ha15_5b.simplices,'sndlit15']

snd_lit16=lithosome(q16,xysnd_knn,400,Heights,'sndlit16')
snd_lit16_b=[x for x in snd_lit16[0] if x[2]>-60 and x[1]>4576020 and x[1]<4578100 and x[0]>425883]
ha16_6b = ConvexHull(snd_lit16_b)
snd_lit16b=[ha16_6b.points,ha16_6b.vertices,ha16_6b.simplices,'sndlit16']

snd_lit17=lithosome(q17,xysnd_knn,500,Heights,'sndlit17')

# List of sands lithosomes

list_sndlithosomesb=[snd_lit1b,snd_lit2b,snd_lit3b,snd_lit4b,snd_lit5b,
                   snd_lit6b,snd_lit7b,snd_lit8b,snd_lit9b,snd_lit10b,
                   snd_lit11b,snd_lit12b,snd_lit13b,snd_lit14b,snd_lit15b,
                   snd_lit16b,snd_lit17]
In [61]:
list_sndlithosomes=[snd_lit1,snd_lit2,snd_lit3,snd_lit4,snd_lit5,
                   snd_lit6,snd_lit7,snd_lit8,snd_lit9,snd_lit10,
                   snd_lit11,snd_lit12,snd_lit13,snd_lit14,snd_lit15,
                   snd_lit16,snd_lit17]

For gravels¶

First, we reshape the data.

In [62]:
len(list_grlithosomes[0])
Out[62]:
5
In [63]:
len(list_sndlithosomes[0])
Out[63]:
5
In [64]:
zip_list_grlithosomes=[list(zip(x[4][0],x[4][1],x[4][2])) for x in list_grlithosomes]

We group the gravel predictions by lithosome.

In [65]:
# data points
list_knn_in_grlit=[[x for x in xyzc_gr50 if ((x[0],x[1],x[2]) in y)] for y in zip_list_grlithosomes]
In [66]:
# confidence values
list_of_c_grlit=[[x[3] for x in y] for y in list_knn_in_grlit]
In [67]:
# heights
z_grlit=[[x[2] for x in y] for y in list_knn_in_grlit]

We work out the average confidence and maximum and minimum heights for each lithosome.

The auxiliar function mean_c2 computes the means.

In [68]:
def mean_c2(lit):
    zlit=list(zip(lit[4][0],lit[4][1],lit[4][2]))
    knlit=[x for x in xyzc_gr50 if ((x[0],x[1],x[2]) in zlit)]
    list_of_c_grlit=[x[3] for x in knlit]
    return np.mean(list_of_c_grlit)
In [69]:
m2=[mean_c2(x) for x in list_grlithosomes]

The following lines perform the same computation as the function mean_c2.

In [70]:
means_cgrlit=[np.mean(x) for x in list_of_c_grlit] # For confidence value
means_zgrlit=[np.mean(x) for x in z_grlit] # For Z coordinate

We check that the same result is obtained.

In [71]:
m2==means_cgrlit
Out[71]:
True

We compute the maximum and minimum height for each lithosome.

In [72]:
max_zgrlit=[max(x) for x in z_grlit]
min_zgrlit=[min(x) for x in z_grlit]

We create the table.

In [73]:
lit_gr_names=['grlit'+str(i+1) for i in range(13)]
In [74]:
d1={'graves lithosomes':lit_gr_names,
    'confidence in mean':means_cgrlit,
    'average height':means_zgrlit,
    'max height':max_zgrlit,
    'min height':min_zgrlit
    }
df = pd.DataFrame(d1)
df.to_csv(DATADIR+'gr_lit_table.csv', index=False, header=True)
In [75]:
df
Out[75]:
graves lithosomes confidence in mean average height max height min height
0 grlit1 50.174238 -72.397590 -40 -100
1 grlit2 51.354275 -51.391417 -40 -60
2 grlit3 51.747555 -55.541150 -45 -70
3 grlit4 49.055077 -56.647059 -50 -65
4 grlit5 50.939566 -45.714286 -35 -50
5 grlit6 49.723211 -28.293173 0 -40
6 grlit7 47.790273 -29.807692 -10 -40
7 grlit8 43.996622 -29.513889 -25 -35
8 grlit9 53.124705 -19.677419 -10 -25
9 grlit10 28.545459 -20.000000 -15 -25
10 grlit11 45.200411 -17.500000 -15 -20
11 grlit12 49.076593 -12.685590 0 -25
12 grlit13 50.428450 -25.000000 0 -40

We can compare the average confidence when taking a grid of size 50 or 200.

In [76]:
gr50_mean # Average confidence for gravel predictions with grid size 50
Out[76]:
49.57830307839629
In [77]:
gr200_mean # Average confidence for gravel predictions with grid size 200
Out[77]:
49.33703116232511

For sands¶

In [78]:
len(list_sndlithosomes)
Out[78]:
17
In [79]:
zip_list_snlithosomes=[list(zip(x[4][0],x[4][1],x[4][2])) for x in list_sndlithosomes]
In [80]:
# data points

list_knn_in_snlit=[[x for x in xyzc_sn50 if ((x[0],x[1],x[2]) in y)] for y in zip_list_snlithosomes]

# We group the sand predictions by lithosome

# Confidence value

list_of_c_snlit=[[x[3] for x in y] for y in list_knn_in_snlit]

# Heights
z_snlit=[[x[2] for x in y] for y in list_knn_in_snlit]

# We work out the average confidence and maximum and minimum heights for each lithosome

means_csnlit=[np.mean(x) for x in list_of_c_snlit]
means_zsnlit=[np.mean(x) for x in z_snlit]

max_zsnlit=[max(x) for x in z_snlit]
min_zsnlit=[min(x) for x in z_snlit]

# Table

lit_sn_names=['snlit'+str(i+1) for i in range(17)]

d2={'sand lithosomes':lit_sn_names,
    'confidence in mean':means_csnlit,
    'cota mean':means_zsnlit,
    'cota max':max_zsnlit,
    'cota min':min_zsnlit
    }


dff = pd.DataFrame(d2)
dff.to_csv(DATADIR+'sn_lit_table.csv', index=False, header=True)
In [81]:
dff
Out[81]:
sand lithosomes confidence in mean cota mean cota max cota min
0 snlit1 50.516362 -83.949275 0 -100
1 snlit2 54.352117 -35.000000 0 -90
2 snlit3 53.833335 -25.512465 0 -90
3 snlit4 52.791327 -39.709507 0 -100
4 snlit5 41.704941 -65.909091 0 -80
5 snlit6 53.765592 -12.724057 0 -60
6 snlit7 54.228871 -24.390244 0 -90
7 snlit8 53.585398 -8.754941 0 -40
8 snlit9 50.545163 -10.555556 -5 -20
9 snlit10 46.070302 -49.139785 -5 -80
10 snlit11 45.906706 -17.142857 0 -40
11 snlit12 54.897272 -10.053763 -5 -55
12 snlit13 42.450885 -40.594059 0 -60
13 snlit14 52.951239 -16.250000 0 -60
14 snlit15 54.124507 -4.559322 0 -10
15 snlit16 53.311568 -9.187117 0 -60
16 snlit17 29.756034 -6.562500 0 -15
In [82]:
sn50_mean # Average confidence for sand predictions with grid size 50
Out[82]:
50.20226163638911
In [83]:
sn200_mean # Average confidence for sand predictions with grid size 200
Out[83]:
50.18378710802425

Drawing the lithosomes¶

In [84]:
# Lithosome data for gravels
data_grlithosomes=[data_lithosome(x[0],x[2],x[3],0,1,'skyblue') for x in list_grlithosomes]

# Lithosome data for sands
data_sndlithosomes=[data_lithosome(x[0],x[2],x[3],0,1,'lemonchiffon') for x in list_sndlithosomesb]
In [85]:
# lithosomes for gravels

dat=data_grlithosomes

fig=go.Figure(data=dat)

fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )))

fig.update_layout( title="3D graves lithosomes Llobregat Delta, Z scale is x 50.",
    scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
                             xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
                             yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
                            ))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_grlithosomes.html',validate=True, auto_open=False)
Out[85]:
'figures/3D_grlithosomes.html'
In [86]:
# lithosomes for sands

dat=data_sndlithosomes

fig=go.Figure(data=dat)

fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )))

fig.update_layout( title="3D sands lithosomes Llobregat Delta, Z scale is x 50.",
    scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
                             xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
                             yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
                            ))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_grlithosomes.html',validate=True, auto_open=False)
Out[86]:
'figures/3D_grlithosomes.html'

Together¶

In [87]:
# lithosomes for gravels and sands

dat=data_grlithosomes+data_sndlithosomes

fig=go.Figure(data=dat)

fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )))

fig.update_layout( title="3D graves and sands lithosomes Llobregat Delta, Z scale is x 50.",
    scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
                             xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
                             yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
                            ))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_lithosomes.html',validate=True, auto_open=False)
Out[87]:
'figures/3D_lithosomes.html'

With basement.

In [88]:
# lithosomes for gravels and sands with basement

dat=data_grlithosomes+data_sndlithosomes

fig=go.Figure(data=dat)

fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
                         x=interpolation_knnbasement50cutted[1],
                         y=interpolation_knnbasement50cutted[2], 
                        opacity = 0.5,
                        #inensity=interpolation_knnbasementcutted[0],
                        colorscale='oryel', 
                        name='basement surface',
                        showscale=False
                        ))

fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )))

fig.update_layout( title="3D graves and sands lithosomes Llobregat Delta with basement surface, Z scale is x 50.",
    scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
                             xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
                             yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
                            ))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_Lithosomes_And_Basement_LRD.html',validate=True, auto_open=False)
Out[88]:
'figures/3D_Lithosomes_And_Basement_LRD.html'

Lithosomes colored by confidence value¶

We modify the function data_lithosome to employ different colors.

In [89]:
def data_lithosome2(vhull,shull,name,alpha,opc,col):
    return go.Mesh3d(x=vhull[:, 0],y=vhull[:, 1], z=vhull[:, 2], 
                        name=name,
                        showlegend=True,
                        colorbar_title=name,
                        color=col, 
                        colorscale='rainbow',
                        i=shull[:, 0], 
                        j=shull[:, 1], 
                        k=shull[:, 2],
                        opacity=opc,
                        alphahull=alpha,
                        showscale=True
                       ) 

Gravels¶

In [90]:
# Data for gravels
data_grlithosomes_conf=[data_lithosome2(x[0],x[2],x[3],0,1,mean_c2(x)) for x in list_grlithosomes]
In [91]:
# the complete gravels lithosome data 

fig=go.Figure(data=data_grlithosomes_conf)

fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )))
fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
                         x=interpolation_knnbasement50cutted[1],
                         y=interpolation_knnbasement50cutted[2], 
                        opacity = 0.5,
                        #inensity=interpolation_knnbasementcutted[0],
                        colorscale='oryel', 
                        name='basement surface',
                        showscale=False
                        ))


fig.update_layout( title="3D graves lithosomes Llobregat Delta colored by confidence, Z scale is x 50.",
    scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
                             xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
                             yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
                            ))
fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_grlithosomes_con.html',validate=True, auto_open=False)
Out[91]:
'figures/3D_grlithosomes_con.html'

Sands¶

In [92]:
len(list_sndlithosomesb)
Out[92]:
17
In [94]:
# Data for sands
data_snlithosomes_conf=[data_lithosome2(list_sndlithosomesb[i][0],
                                    list_sndlithosomesb[i][2],
                                    list_sndlithosomesb[i][3],
                                    0,1,mean_c2(list_sndlithosomes[i])) for i in range(17)]
In [95]:
fig=go.Figure(data=data_snlithosomes_conf)

fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )))

fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
                         x=interpolation_knnbasement50cutted[1],
                         y=interpolation_knnbasement50cutted[2], 
                        opacity = 0.5,
                        #inensity=interpolation_knnbasementcutted[0],
                        colorscale='oryel', 
                        name='basement surface',
                        showscale=False
                        ))


fig.update_layout( title="3D sand lithosomes Llobregat Delta colored by confidence, Z scale is x 50.",
    scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
                             xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
                             yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
                            ))

fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_snlithosomes_con.html',validate=True, auto_open=False)
Out[95]:
'figures/3D_snlithosomes_con.html'

Together with basement.

In [96]:
# Data for gravels with opacity 0.3
data_grlithosomes_conf05=[data_lithosome2(x[0],x[2],x[3],0,0.5,mean_c2(x)) for x in list_grlithosomes]
In [97]:
dat=data=data_snlithosomes_conf+data_grlithosomes_conf

fig=go.Figure(data=dat)
fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
                         x=interpolation_knnbasement50cutted[1],
                         y=interpolation_knnbasement50cutted[2], 
                        opacity = 0.5,
                        #inensity=interpolation_knnbasementcutted[0],
                        colorscale='oryel', 
                        name='basement surface',
                        showscale=False
                        ))

fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )))

fig.update_layout( title="3D gravels and sand lithosomes Llobregat Delta colored by confidence, Z scale is x 50.",
    scene=dict(aspectratio=dict(x=2, y=2, z=0.5),
                             xaxis = dict(range=[bound[0]-2000,bound[1]+2000],),
                             yaxis = dict(range=[bound[2]-2000,bound[3]+2000])
                            ))

fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_lithosomes_con.html',validate=True, auto_open=False)
Out[97]:
'figures/3D_lithosomes_con.html'

Comparing gravels¶

We will now display two figures side by side. One with the lithosomes in a single color and another with lithosomes colored by the average confidence.

We need to import subplots from plotly.

In [98]:
from plotly.subplots import make_subplots

Without basement.

In [99]:
fig = go.Figure().set_subplots(1, 2, horizontal_spacing=0.05,
                               specs=[[{"type": "scatter3d"}, {"type": "scatter3d"}]],
                               subplot_titles=('graves lithosomes','graves lithosomes colored by confidence')
                              )

for i in range(len(data_grlithosomes)):
     fig.add_trace(data_grlithosomes[i],row=1, col=1)
        

#fig.add_trace(marks_data[0],row=1, col=1)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           legendgroup = '1',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )), row=1, col=1)

#fig.update_layout(legend=dict(x=0,y=1))
    
for i in range(len(data_grlithosomes_conf)):
    fig.add_trace(data_grlithosomes_conf[i],row=1, col=2)


fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )), row=1, col=2)

fig.update_traces(showlegend=True)
fig.update_layout( title="3D graves lithosomes comparation Llobregat Delta, Z scale is x 50.",
                  scene1=dict(aspectratio=dict(x=1, y=1, z=0.5)
                            #aspectmode='auto',
                            #aspectratio=dict(x=2, y=2, z=1)
                            ),
                   scene2=dict(aspectratio=dict(x=1, y=1, z=0.5)
                            #aspectmode='auto',
                            #aspectratio=dict(x=2, y=2, z=1)
                            ),
                  legend=dict(x=1,y=0),
                  height=600, showlegend=True)

#fig.update_layout(scene2_aspectmode='manual',
 #                 scene2_aspectratio=dict(x=1, y=1, z=0.5),
  #                legend=dict(x=1,y=0)
   #              )
    

fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_grlithosomes_comparation_nobasement.html',validate=True, auto_open=False)
Out[99]:
'figures/3D_grlithosomes_comparation_nobasement.html'

With basement.

In [100]:
fig = go.Figure().set_subplots(1, 2, horizontal_spacing=0.05,
                               specs=[[{"type": "scatter3d"}, {"type": "scatter3d"}]],
                               subplot_titles=('graves lithosomes','graves lithosomes colored by confidence')
                              )

for i in range(len(data_grlithosomes)):
     fig.add_trace(data_grlithosomes[i],row=1, col=1)
        
fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
                         x=interpolation_knnbasement50cutted[1],
                         y=interpolation_knnbasement50cutted[2], 
                        opacity = 0.5,
                        #inensity=interpolation_knnbasementcutted[0],
                        colorscale='oryel', 
                        name='basement surface',
                        showscale=False
                        ), row=1, col=1)
#fig.add_trace(marks_data[0],row=1, col=1)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           legendgroup = '1',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )), row=1, col=1)

#fig.update_layout(legend=dict(x=0,y=1))
    
for i in range(len(data_grlithosomes_conf)):
    fig.add_trace(data_grlithosomes_conf[i],row=1, col=2)
#fig.add_trace(marks_data[0],row=1, col=2)
fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
                         x=interpolation_knnbasement50cutted[1],
                         y=interpolation_knnbasement50cutted[2], 
                        opacity = 0.5,
                        #inensity=interpolation_knnbasementcutted[0],
                        colorscale='oryel', 
                        name='basement surface',
                        showscale=False
                        ), row=1, col=2)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )), row=1, col=2)

fig.update_traces(showlegend=True)
fig.update_layout( title="3D graves lithosomes comparation Llobregat Delta, Z scale is x 50.",
                  scene1=dict(aspectratio=dict(x=1, y=1, z=0.5)
                            #aspectmode='auto',
                            #aspectratio=dict(x=2, y=2, z=1)
                            ),
                   scene2=dict(aspectratio=dict(x=1, y=1, z=0.5)
                            #aspectmode='auto',
                            #aspectratio=dict(x=2, y=2, z=1)
                            ),
                  legend=dict(x=1,y=0),
                  height=600, showlegend=True)

#fig.update_layout(scene2_aspectmode='manual',
 #                 scene2_aspectratio=dict(x=1, y=1, z=0.5),
  #                legend=dict(x=1,y=0)
   #              )
    

fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_Gravel_Lithosomes_Confidence.html',validate=True, auto_open=False)
Out[100]:
'figures/3D_Gravel_Lithosomes_Confidence.html'

Comparing sands¶

Without basement.

In [101]:
fig = go.Figure().set_subplots(1, 2, horizontal_spacing=0.05,
                               specs=[[{"type": "scatter3d"}, {"type": "scatter3d"}]],
                               subplot_titles=('sands lithosomes','sands lithosomes colored by confidence')
                              )

for i in range(len(data_sndlithosomes)):
     fig.add_trace(data_sndlithosomes[i],row=1, col=1)
#fig.add_trace(marks_data[0],row=1, col=1)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )), row=1, col=1)


for i in range(len(data_snlithosomes_conf)):
    fig.add_trace(data_snlithosomes_conf[i],row=1, col=2)
#fig.add_trace(marks_data[0],row=1, col=2)
fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )), row=1, col=2)

#fig.update_layout(legend=dict(
 #   yanchor="top",
  #  y=0.99,
   # xanchor="right",
    #x=0.01))

fig.update_layout( title="3D sands lithosomes comparation Llobregat Delta, Z scale is x 50.",
                  scene1=dict(aspectratio=dict(x=1, y=1, z=0.5)
                            #aspectmode='auto',
                            #aspectratio=dict(x=2, y=2, z=1)
                            ),
                   scene2=dict(aspectratio=dict(x=1, y=1, z=0.5)
                            #aspectmode='auto',
                            #aspectratio=dict(x=2, y=2, z=1)
                            ),
                  height=600, showlegend=True)
    
                  

fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_snlithosomes_comparation_nobasement.html',validate=True, auto_open=False)
Out[101]:
'figures/3D_snlithosomes_comparation_nobasement.html'

With basement.

In [102]:
fig = go.Figure().set_subplots(1, 2, horizontal_spacing=0.05,
                               specs=[[{"type": "scatter3d"}, {"type": "scatter3d"}]],
                               subplot_titles=('sands lithosomes','sands lithosomes colored by confidence')
                              )

for i in range(len(data_sndlithosomes)):
     fig.add_trace(data_sndlithosomes[i],row=1, col=1)
#fig.add_trace(marks_data[0],row=1, col=1)

fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
                         x=interpolation_knnbasement50cutted[1],
                         y=interpolation_knnbasement50cutted[2], 
                        opacity = 0.5,
                        #inensity=interpolation_knnbasementcutted[0],
                        colorscale='oryel', 
                        name='basement surface',
                        showscale=False
                        ), row=1, col=1)

fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )), row=1, col=1)


for i in range(len(data_snlithosomes_conf)):
    fig.add_trace(data_snlithosomes_conf[i],row=1, col=2)
#fig.add_trace(marks_data[0],row=1, col=2)

fig.add_trace(go.Surface(z=interpolation_knnbasement50cutted[0],
                         x=interpolation_knnbasement50cutted[1],
                         y=interpolation_knnbasement50cutted[2], 
                        opacity = 0.5,
                        #inensity=interpolation_knnbasementcutted[0],
                        colorscale='oryel', 
                        name='basement surface',
                        showscale=False
                        ), row=1, col=2)

fig.add_trace(go.Scatter3d(x=xyzcontour[0],y=xyzcontour[1],z=xyzcontour[2], mode="lines",
                           line_width=5,
                           name='Delta Contour',
                           marker = dict(
                               size = 4,
                               color = 'black'
                               )), row=1, col=2)

#fig.update_layout(legend=dict(
 #   yanchor="top",
  #  y=0.99,
   # xanchor="right",
    #x=0.01))

fig.update_layout( title="3D sands lithosomes comparation Llobregat Delta, Z scale is x 50.",
                  scene1=dict(aspectratio=dict(x=1, y=1, z=0.5)
                            #aspectmode='auto',
                            #aspectratio=dict(x=2, y=2, z=1)
                            ),
                   scene2=dict(aspectratio=dict(x=1, y=1, z=0.5)
                            #aspectmode='auto',
                            #aspectratio=dict(x=2, y=2, z=1)
                            ),
                  height=600, showlegend=True)
    
                  

fig.show()
go_offline.plot(fig,filename=FIGURESDIR+'3D_Sand_Lithosomes_Confidence.html',validate=True, auto_open=False)
Out[102]:
'figures/3D_Sand_Lithosomes_Confidence.html'
In [ ]: